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Abstract. In this paper we investigate the effect of stochasticity in the spatial and temporal 
distribution of supernova remnants on the spectrum and chemical composition of cosmic 
rays observed at Earth. The calculations are carried out for different choices of the diffusion 
coefficient D{E) experienced by cosmic rays during propagation in the Galaxy. In particular, 
at high energies we assume that D{E) oc , with 5 = 1/3 and 6 = 0.6 being the reference 
scenarios. The large scale distribution of supernova remnants in the Galaxy is modeled 
following the distribution of pulsars, with and without accounting for the spiral structure of 
the Galaxy. We find that the stochastic fiuctuations induced by the spatial and temporal 
distribution of supernovae, together with the effect of spallation of nuclei, lead to mild but 
sensible violations of the simple, leaky-box-inspired rule that the spectrum observed at Earth 
is N{E) oc E^"' with a = ^ + 6, where 7 is the slope of the cosmic ray injection spectrum 
at the sources. Spallation of nuclei, even with the small rates appropriate for He, may 
account for small differences in spectral slopes between different nuclei, possibly providing 
an explanation for the recent CREAM observations. For 6 = 1/3 we find that the slope of the 
proton and helium spectra are ~ 2.67 and ~ 2.6 respectively (with fiuctuations depending 
on the realization of source distribution) at energies around ~ 1 TeV (to be compared with 
the measured values of 2.66 it 0.02 and 2.58 it 0.02). For 6 = 0.6 the hardening of the He 
spectra is not observed. The stochastic effects discussed above cannot be found in ordinary 
propagation calculations, such as GALPROP, where these effects and the point like nature 
of the sources are not taken into account. We also comment on the effect of time dependence 
of the escape of cosmic rays from supernova remnants, and of a possible clustering of the 
sources in superbubbles. In a second paper we will discuss the implications of these different 
scenarios for the anisotropy of cosmic rays. 
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1 Introduction 

Any investigation of the origin of cosmic rays (CRs) has to address some basic questions: 
1) where and how are CRs produced/accelerated? 2) what is their chemical composition at 
the sources? 3) How do they propagate to us? 4) What is the anisotropy that results from 
source distribution (in space and time) and from propagation? 

After the pioneering proposal of Baade and Zwicky [1], the idea that the bulk of Galactic 
CRs may be produced in supernova remnants (SNRs) has become increasingly more popular. 
The main reason for this success is that on energetic grounds these astrophysical sources are 
the most reliable possibility in that an efficiency of ~ 10% in the form of accelerated particles 
seems to be able to provide a qualitatively good description of the CR fluxes observed at 
Earth. The second major reason for the success has come several decades after the original 
idea was put forward, and relies on the proposal of diffusive shock acceleration (DSA) [2, 3, 4] 
as a conversion mechanism from kinetic energy of the expanding supernova blast wave into 
accelerated particles. This mechanism is seen at work in shocks in the solar system and 
is based on relatively simple principles. These two ingredients, together with a bulk of 
observations, showing likely hints of efficient CR acceleration in SNRs, have determined the 
upgrading of the SNR idea to the rank of a standard paradigm. Yet, the paradigm has not 
been proven right beyond doubt so far, as discussed recently in Ref. [5]. 

The problems with confirming the paradigm originate primarily from addressing all the 
four questions reported above self-consistently: on one hand the spectra of CRs accelerated 
in SNRs, according with test-particle DSA are very close to power laws N{E) cx E~''' with 
slope 7 ~ 2. On the other hand, the spectra measured at Earth require a diffusion coefficient 
D{E) DC with 5 ~ 0.7, so that the equilibrium spectrum is n[E) oc N[E)/D{E) cc 
^-7-(5 _ ^-2.7 _ However the strong diffusion that follows from this line of thought seems to 
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lead to exceedingly large anisotropy [6]. Non linear effects in DSA make this problem even 
more severe in that spectra are predicted to be concave and lead to require an even stronger 
energy dependence of the diffusion coefficient at high energies, as found in Ref. [7], where, 
however, the problem of anisotropy was not discussed. 

It has recently been proposed that the non-linear theory of DSA (hereafter NLDSA) 
may lead to somewhat steeper spectra if the velocity of the scattering centers is taken into 
account in solving the transport equation for CRs at SNR shocks [8, 9, 10]. In this case 
the spectra of accelerated particles may be less concave, and have an approximate power 
law shape with slope ~ 2.1 — 2.2, thereby implying 6 ~ 0.5 — 0.6. The steeper slope of the 
injection spectrum follows from both the finite velocity of the scattering waves and from 
the convolution over time of the acceleration history of SNRs throughout their evolution in 
the interstellar medium (ISM). It is noteworthy, however, that this mechanism introduces a 
very disappointing element in that the result depends rather strongly on poorly understood 
details such as the wave helicity, as well as the reflection and transmission of waves at the 
shock surface (in principle the effect considered here could even lead to harder spectra). 

The most serious limitation to the identification of SNRs as the sources of Galactic CRs 
is the lack of a firm detection of gamma rays that can unequivocally be interpreted as due 
to production and decay of neutral pions: in the few cases in which this association appears 
to be easier, the SNR is close to a molecular cloud (denser target for pp interactions) and 
typically of old age. In such circumstances the spectrum of accelerated particles is unlikely 
to be representative of SNRs at their best in terms of CR production. In these few cases, 
observations show a rather steep spectrum of accelerated particles, quite unlike the ones 
predicted by NLDSA. Possible exceptions are RX J1713. 7-3946 and Tycho, although for the 
first remnant other problematic aspects may be found [11, 12]. Moreover, a recent analysis of 
the data on this remnant by the Fermi collaboration [13] suggests that in fact the spectrum 
of gamma rays is very flat thereby leading to conclude that the emission is likely of leptonic 
origin. On the other hand the recent gamma ray detection of the Tycho supernova remnant 
by the VERITAS Cherenkov telescope [14] and by the Fermi telescope [15] appear to be best 
fit by a model in which particle acceleration occurs efficiently [16]. 

As discussed in [9, 10], a qualitatively good fit to the spectra of CRs (including nuclei 
heavier than hydrogen) can be achieved within the context of NLDSA with finite speed of 
scattering centers (implying 6 ~ 0.55) but some spectral features that have been recently 
found, the so-called discrepant hardenings [17] observed by CREAM, are not explained as 
yet within the context of the SNR paradigm (for a different, environment related explanation, 
see [18]). 

From the point of view of propagation in the Galaxy, the standard calculations carried 
out with GALPROP [19], or similar propagation codes, suggest that one can obtain a good 
global fit also with reacceleration models, namely with diffusion coefficient D{E) oc E^^'^ and 
second order Fermi acceleration: the latter is relevant only at low energies, where it allows 
one to better reproduce the B/C ratio. This scenario would imply a high energy injection 
spectrum N{E) oc E^^''^, very challenging for NLDSA in SNRs, but certainly preferable from 
the point of view of anisotropy (see Paper II [20]). 

In the present work we use the Green function formalism to describe the propagation 
of CRs to Earth if the sources are SNRs. These are modeled as discrete sources in space 
and time, with a spatial distribution and a rate deduced from observations. The effect of 
discreteness on the spectrum and chemical composition of CRs observed at Earth is the goal 
of this first paper. In a second paper [20] we will discuss the more delicate issue of anisotropy. 
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We concentrate here on 1) the investigation of the spectra and chemical fluctuations induced 
by the discrete nature of the sources; 2) the chemical composition around the knee; 3) the 
spectral differences between protons and heavier elements; 4) the effect of extragalactic CRs 
on the end of the Galactic CR spectrum. 

Previous attempts at taking into account the random nature of the sources on the CR 
spectrum observed at Earth were presented in [6, 21]. The first paper mainly focuses on the 
issue of anisotropy and will be discussed in more detail in the forthcoming Paper II. 

In [21], the authors use a series expansion of the propagation equation to show that the 
stochasticity of source distribution induces large temporal fluctuations in the spectrum of 
primary CR nuclei. Assuming a diffusion coefficient D{E) oc E'^'^, they find fluctuations by 
20% on average, but up to 100%, which leads them to question the goodness of the B/C ratio 
as an indicator of CR propagation parameters in the Galaxy. However, a diffusion coefficient 
scaling as E^'^ implies a level of anisotropy well in excess of the observations, as discussed in 
[6] and in Paper II. If a weaker energy dependence of the diffusion coefficient is considered, 
the fluctuations found by [21] are much reduced, weakening also their conclusion. 

The present paper is organized as follows: in § 2 we introduce the relevant Green 
functions; in § 3 we derive some basic results for the simple case of a homogeneous distribution 
of sources in the disc of the Galaxy and discuss the limitations of the diffusion formalism for 
individual sources in § 3.1. In § 4 we derive simple estimates of the effect of fluctuations on 
the spectrum for a homogeneous but discrete distribution of sources in the disc. Our results 
on spectra and chemical composition for a stochastic distribution of the sources in the disc 
are presented in § 5), where the spiral structure of the Galaxy is first ignored (§ 5.1) and 
then taken into account (§ 5.2). Additional scenarios are discussed in § 6, while we discuss 
the transition region between galactic and extragalactic CRs in § 7. We conclude in § 8. 

2 Relevant Green functions 

The diffusive transport of cosmic rays from a point source located at a position rg = {xsiUs, ^s) 
and injecting a spectrum N{E) at a time tg can be written as: 



where nk{E,f,t) is the density of particles of type k (nuclei) with energy E at the location 
r and time t, Dk{E) is the diffusion coefficient assumed to be spatially constant and T^j^{E) 
is the rate of spallation of nuclei of type k to lead to lighter nuclei. Each source is assumed 
to produce nuclei of H (fc = 1), He (fc = 2), CNO (A: = 3), Mg-Al-Si {k = 4) and Fe (A; = 5). 
The cutoff energy of the proton component, E^^^^ is taken in such a way as to achieve the 
best fit to the all-particle CR spectrum. For heavier nuclei the maximum energy is chosen 
to be Z times larger than for H, where Z = 2 for He, Z = 7 for CNO, Z = 13 for Mg-Al-Si 
and Z = 26 for Fe nuclei. All injection spectra are therefore in the form: 



where the proportionality constant for each of the nuclear species is determined a posteriori 
in order to fit the normalization of the spectra of nuclei of type k as measured at Earth by 
CREAM at the energy of 1 TeV [17, 22]. It is worth recalling that the scaling of the maximum 
energy with the charge Z of the nucleus relies upon the assumption of full ionization of the 



dnk{E,f,t) 
dt 



V[Dk{E)Vnk{E,f,t)] 



rf (E)nfc(E, f, t) + Nk{E)6{t - ts)d'^{r 



rs), (2.1) 




(2.2) 
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nucleus. In fact, this condition is achieved only during the acceleration process and the 
actual maximum energy of heavy nuclei might be less than Z times the maximum energy of 
protons, because of the difficulty in achieving photo-ionization of atoms through extraction 
of electrons in the inner atomic shells [23]. 

The diffusion region is assumed to be a cylinder of infinite radius and half height H. 
Inside this cylinder the diffusion coefficient is constant. The fact that the radius is assumed 
to be potentially infinite is not a limitation to the calculation because we will see that the 
source distributions that will be adopted are concentrated in a region of 4 — 10 kpc radius 
(the Galactic disc) and for practical purposes the size of the diffusion region can be taken 
to be much larger than the radius of the disc. The escape of CRs from this model Galaxy 
occurs through the upper and lower boundaries. The escape is modeled by assuming that 
71/; vanishes at z = ±i?, and the escape fiux through the surfaces z = zizH is described by 

In this paper we concentrate our attention on particles with energy above ~ 1 TeV, 
therefore in Eq. 2.1 we ignored terms of advection and terms of possible second order reac- 
celeration which are both potentially important at much lower energies. We also neglect 
the contribution to the flux of nuclei deriving from secondary nuclei produced in spallation 
events. From the point of view of diffusive propagation the calculation presented here is not 
too different from GALPROP or similar propagation codes. The only two relevant differ- 
ences are in the more phenomenological way in which spallation is described and in the fact 
that here we do not impose the free escape boundary condition at some finite radius along 
the lateral sides of the cylinder. The latter, as discussed above, is not a relevant limitation 
although it might lead to some corrections for large values of H. The former is not important 
in the present context since we will limit ourselves to primary nuclei, while it might represent 
a limitation if we were to describe rare nuclei, such as B and ^^Be or similar secondary 
products. 

The spallation rate is defined in terms of the gas density in the diffusion region and of the 
cross section for spallation: T'JF'''{E) = UgasCagp^k (the velocity of nuclei has been assumed to 
equal the speed of light c since all our calculations will only concern ultra-relativistic nuclei). 
The cross section for spallation of a nucleus k (mass A^) can be written as [24]: 

(Tsp,k{E) = asp{E)A{^-^''^ mb, (2.3) 

where 

asp{E) = 50.44 - 7.931og(£;ey) + 0.61 [\og{E^v)? (2.4) 



and 



Psp{E) = 0.97 - 0.022 \og{E^v) (2.5) 



where Ef^y is the energy of the nucleus of type k in units of eV. The gas density is 
non-trivial to define in a homogeneous model as this is. We assume here that the gas density 
in the disc of the Galaxy (with half-thickness h) is n^igc = 3 cm~^ while it vanishes in the 
halo. In this case the mean density experienced by CRs during propagation is approximately 

ngas ~ ndisc{h/H). 

The free Green function (namely without boundary conditions) of Eq. 2.1 can be written 

as 



Qi''^\r,t]fs,ts) = -^^^^exp [-rf (^)r] exp 
[47rDfcr] ' 



(r 



ADkT 



(2.6) 
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where r = t — tg- The Green function that satisfies the correct boundary condition at z = zizH 
can be obtained by using the image charge method and can be written as follows: 



[4TTDkTf 



Nk{E) 



3/2 



exp [-T'IJ'{E)t] exp 



(x - Xsf + {y- Vs) 



21 



X 




E (-ir^xp 



(2.7) 



n=—oo 



where z'^ = (— 1)"^^ + 2nH are the z coordinates of the image sources. It is easy to check 
that Qk{x, y,z = m, t; Xs,ys, Zs,ts) = by simply expanding the sum term. In what follows 
the Earth will be located at {x, y, z) = (i?©, 0, 0) with Rq = 8.5 kpc the distance of the Sun 
from the center of the Galaxy. It is however useful to keep Eq. 2.7 in its most general form 
because anisotropics are related to the spatial derivatives of the Green function calculated 
at the detection location. 

It is important to realize that the Green function formalism outlined here allows us to 
take into account a completely arbitrary spatial distribution of the sources in the Galaxy, 
and also an arbitrary temporal evolution of the injection of cosmic rays from each source. 

Two instances related to the specific case of SNRs as sources of CRs may clarify the 
importance of this point. SNRs produce CRs in a rather complicated and not entirely clear 
way: while the spectrum of CRs accelerated at the external shock is well defined and can 
be calculated using the theory of NLDSA, the spectrum of particles escaping a remnant to 
become CRs is much more troublesome in that it depends on how and when CRs escape from 
the acceleration region. Most accelerated particles in a SNR shock are advected downstream 
of the shock and are trapped in the expanding shell for the whole duration of the Sedov- Taylor 
phase, thereby losing energy adiabatically. These particles become CRs in the ISM only after 
the SNR ends its evolution and the shock dies out. On the other hand, during the Sedov- 
Taylor phase and in presence of magnetic field amplification and damping, the maximum 
energy of accelerated particles decreases with time. At any given time the particles with 
the instantaneous maximum energy may escape the system from upstream. The spectrum 
of these particles at any given time t is close to a delta function 6{E — Emaxit)), where 
Emax{t) is the maximum energy reached at that time. The convolution over time of these 
peaked spectra leads to a power law injection spectrum [25] which is not directly related to 
the spectrum of particles at the shock. The global injected spectrum from an individual SNR 
is likely to be the result of the superposition of the spectrum of particles escaping the SNR 
from upstream and of those that escape after the end of the Sedov- Taylor phase. The whole 
process lasts ~ (3 — 30) x 10^ years, with the low energy particles escaping at later times 
than the higher energy ones. During the ejecta dominated phase that precedes the adiabatic 
phase, the maximum energy of accelerated particles grows and escape is not effective. 

In our calculations we consider two different models of injection: 1) burst injection (all 
particles from a SNR are injected at the same time with a spectrum N(E)); 2) continuous 
injection in which at any time the injected spectrum is oc 6[E — -Emax(i)]- This second choice 
should qualitatively reflect the statement that escape of particles from the accelerator takes 
place in such a way that higher energy particles escape earlier in the evolution, while lower 
energy particles leave later. 
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In the first scenario, the Green function in Eq. 2.7 is already the solution that we are 
seeking. The spectrum injected by an individual SNR can be written in the form 



NkiE) 



(7fe-2)r/fcefcm 1 / E \ f E 



1 



"7fe+2 



^0 / \ ^max,k , 



where rjk is the fraction of the kinetic energy of the blast wave, ekin^ that goes into accelerated 
particles of type k. The reference energy £o,fc is taken to be 1 GeV for protons while for 
heavier nuclei its numerical value is not really important since we simply rescale the injected 
protons spectrum in order to fit the spectra observed at Earth. In Eq. 2.8, Emax,k is the 
maximum energy of particles of type k. In the expression above and in what follows we always 
assume that the injection spectrum is steeper than E~'^, since flatter spectra would result in 
unreasonable choices of the diffusion coefficient, namely in exceedingly large anisotropy and 
even in the breaking of the regime of diffusive propagation (see below). 

The second model of injection that we consider is introduced, as we already mentioned, 
in order to take into account, at least in a qualitative way, the fact that escape of particles 
from the accelerator is expected to occur in a differential way, with higher energy particles 
escaping first. As discussed in [8, 25], the escape starts at the beginning of the Sedov- 
Taylor (ST) phase, while during the ejecta dominated (ED) phase of the SNR evolution the 
particles are confined in the expanding shell, the maximum energy increases with time and 
the particles lose energy adiabatically. During the ST phase, in the presence of magnetic field 
amplification, the maximum energy decreases after reaching its maximum value between the 
end of the ED phase and the beginning of the ST phase. 

The ST phase starts when the mass swept up by the supernova shock equals the mass 
of the ejecta, Mej: 

Me,- = -t^pR\t ^ RsT = ( ^ ) = 6-6 X lO^^ ( ) c"^' (2-9) 

where Rst is the radius of the shell at the beginning of the ST phase, M^j^q is the mass of 
the ejecta in units of solar masses and ni is the gas density in the ISM in units of 1 particle 
per cubic cm. The kinetic energy of the explosion is related to the mass of the ejecta as 

Cfcin = ^MejUsT (2.10) 

which leads to 

where ust is the velocity of the expanding shell at the beginning of the ST phase. It is 
interesting to notice that 

ekin = ^MejulT = (^p4t) {^-^Rst) ^st ■ (2.12) 

This equation, together with Eq. 2.11, gives an expression for the time Tst at which the ST 
phase approximately starts: 

TsT = \^^70yr^^. (2.13) 

3 UST fl( n}' 
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During the ST phase the radius and velocity of the sheh change as 



Rsh{t) = RsT 



2/5 



Ush{t) 



tUsT 



TsT 



-3/5 



(2.14) 



We assume that escape of CRs from the SNR occurs from upstream of the shock, 
therefore the spectrum of escaping particles at any given time is very much peaked around 
Emax{t) (we omit here the index k that identifies the type of nucleus) and for simplicity 
we assume Q{E,t) = K5 [E — Emax{t)\- In the absence of a fundamental description of the 
way Emax changes in time (it depends on the details of magnetic field amplification and 
damping) we assume, for t > Tst, Emax{t) = EM{t/TsT)~°^ (a > 0), with Em the maximum 
energy reached at the time Ts- The value of a is chosen in such a way as to achieve that the 
maximum energy at the time when the SNR dies out, tsnr, is Emax{TSNR) = Eq where we 
take Eq = 1 GeV (for protons). The normalization constant K is calculated by integrating 
Q{E) over E: 



dEQ{E,t)E = r]{t) 



(47ri2,\) ^ K(t) = vit) 



1 



Tst Em \Tst 



t 



(2.15) 

For reasons that will be clear in the following, we also assume that r/(t) = r]Q{t/TsT)^ , 
with P >Q. It follows that the density of CRs from a SNR at location is 



n(£,f,t) 



Min[TsT+rsNR,t-ts\ 



dt*7]{t* 



1 



t* 



E — Em 



Tst \5J Em \Tst 



gk{r,t;fs,t*) 



a-l 



(2.16) 



The variable r in Eq. 2.7 is clearly r = t — t* — ts, and the (5-function in Eq. 2.16 selects 
the time t* = Tst{Em / E)^/'^ . One can then write: 



n{E, f, t) = ?7o ek 
1 



aEl 



E 



\ E 



[ATrD{E)T*f^ 



exp [-rf (i?)r*] exp 

+ 00 



{x - Xsf + {y- Vs) 



21 



/ \2 



- 4) 



(2.17) 



where 



Em 



E 



T*=t-ts- Tst 
From Eq. 2.16 it is clear that the solution is non-zero only for 

Tsnr t — ts 



1<(^Y <Min 



1 + 



Tst ' Tst 
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The first inequality is satisfied by definition of Em- The second inequality reads differently 
for old and young SNRs. For old SNRs, namely the ones that have completed their evo- 
lution before the observation time t, the condition reads E > Em ^1 + ^Tg^ ) ' which is 
again satisfied by definition for the energies we are interested in, since we required that 
Em ^ ^TsT ) ~ interesting is the case in which (t — ts)/TsT < 1 + tsnr/Tst, 

namely the case of recent SNRs. In this case the spectrum in Eq. 2.17 is non-zero only for 

For instance, for tsnr = 10^ years, Em = 3 x 10^ GeV and Tst ~ 1000 years one has a ^ 3.2 
in order to satisfy Em{tsnr) = Eq =lGeV. Therefore from a supernova that went off 20,000 
years ago, such as Vela, we can only receive particles with energy E > 200 GeV. It is worth 
noticing that this fact has nothing to do with the propagation time, which adds to limiting 
the minimum energy of the particles that can be detected at Earth. The limit discussed here 
is due instead to the fact that low energy particles could not leave the source even at relatively 
long times after the beginning of the ST phase. This effect is unlikely to be important for 
the total spectrum of CRs observed at Earth, which is dominated by the superposition of 
numerous distant SNRs, but it can potentially be important for the anisotropy signal as we 
discuss in Paper II. 

One last point that is worth noticing about Eq. 2.17 is related to the injection spectrum 
that can be obtained in the approach illustrated above: one can see that if the acceleration 
efficiency is constant in time (^5 = 0) the spectrum is always oc E~'^ and is totally independent 
on the spectrum of particles accelerated at the shock. The power law shape of the spectrum 
is uniquely related to the evolution in time of the SNR and not to the acceleration process 
at work at the shock. This issue, discussed at length in [8], is hard to evade and in fact 
all physical effects that can be added to this approximate picture almost invariably lead 
to spectra even flatter than E~'^. Therefore a crucial assumption that is required in order 
to have injected spectra steeper than E~'^ is that /3 > 0. Physically, this implies that the 
acceleration efficiency is required to increase with time while the SNR gets older. At present 
it has not been possible to find any physical justification for this requirement. 



3 Simple benchmark results on spectra 

In this section we provide a benchmark calculation of the spectrum of nuclei expected from 
discrete sources (SNRs) in the disc of the Galaxy. The full calculation, taking into account 
the spatial and temporal distribution of SNRs in the Galaxy, will be illustrated in § 5. Here 
we simply present an analytical calculation that shows the main relevant scaling relations. 
For this purpose we consider the disc as infinitely thin in the z direction (the sources are all 
located in a plane with Zs = 0), and having a radius Rd- The flux of CRs is calculated at the 
center of the disc. In a homogeneous model as this is, this does not induce a large error. The 
spectrum of protons (no spallation) for bursting sources exploding with a rate TZ is easily 
written as 



ncR{E) 



dr 



R. 



27rr 
dr- 



N{E)n 



T^Rd [AttD{E)t 



3/2 



exp 



AD{E)t 



E (-l)"exp 



{2nHf 
'AD{E)t 

(3.1 
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Carrying out the integration on r first and then on r, one easily obtains 



ncR{E) 



N{E)n 
27rD{E)Rd 



E(- 



1 + 



/2nH 



/2nH 



(3.2) 



One can check that the sum over n equals ~ H/R^ for H <^ so that 



ncR{E) 



N{E)n H _ N{E)n 



2^Rl 



D{E) 2H7TRf, D{E)' 



(3.3) 



The first expression clarifies that the flux of CRs scales with H/D{E), a result that remains 
valid even in more complex versions of the calculation and in fact found also in propagation 
codes such as GALPROP [19] and DRAGON [26]. The second expression in Eq. 3.3 clarifies 
that the observed density of CRs is simply equal to the total injection rate N{E) TZ divided 
by the volume of the Galaxy 2H7rR'^, and multiplied by the escape time H"^ / D{E). 

It is immediately apparent that the diffusion coefficient plays a crucial role in this type 
of calculations. Since in this paper we do not calculate the B/C ratio or the abundance 
of radioactive isotopes, we choose to borrow the normalization of D{E) from the results of 
propagation codes such as GALPROP and DRAGON. They agree on the conclusion that, if 
the diffusion coefficient is chosen to be in the form 



D{E) = 10^^ L>: 



28 



^ \ 2-1 



(3.4) 



for rigidity R > 2> GV, the secondary to primary ratios and the abundances of unstable 
isotopes can be best fit by choosing D2s/Hkpc = 1-33 for 5 = 1/3 and D2^/Hkpc = 0.55 for 
5 = 0.6, where H^pc is the height of the halo in units of kpc. 

Let us now discuss the case of nuclei, which is somewhat more interesting due to the 
effect of spallation. The density of CR nuclei of type k in our simple model can be written 
as ^ 

nk{E) = I dr [ dr^^ ^ ^J^^^l^.-xi-, ^^p 







exp 



ADk{E)T 



T^Rj [47TDkiE)r] 

+ 00 

E (-i)"exp 



'sp,k 



{2nHf 
'^Dk{E)T 



(3.5) 



where we defined Tsp.k = ^/^sp,k- The integrals over r and r are both carried out analytically 
in this order, leading to: 

Nk{E)n 



nk{E) 



27r Dk{E)Rj 



Dk{E)Tsp,kX 



+ 00 



exp 



An 



2 '^esc,k 
'^sp,k 



1/2 



exp 



Tesc,k 
"Tspfk 



1/2 



r2 

4^ +^2 



1/2- 



(3.6) 



where Tesc,k{E) = H"^ / Dk{E) is the escape time of nuclei of type k. The sum in Eq. 3.6 
tends to ■\/Tesc,k/rsp,k for Tesc,k/Tsp,k <■ 1 and to 1 for Tesc,k/Tsp,k > 1. It follows that when 
spallation is negligible 



nk{E) 



Nk{E)n H Tesck 

27rRl Dk{Ey Tsp,k 



<C 1 



(3.7) 
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which is the same solution as for protons (Eq. 3.3), showing an energy dependence nk{E) cx 
^-7-<5_ ^Yie opposite hmit (spallation dominant over escape) one has to be careful in 
defining the effective density experienced by the propagating CRs, which now becomes energy 
dependent and given by 

For the CR spectrum one obtains: 

/px Nk{E)TZ Tesc,k ^ . „x 

^k[E) = OLJ d2 V'^sp,kTesc,k , > 1, (3-9) 

which leads to ni.{E) oc E~'^ , due to the fact that with Ugas given by Eq. 3.8, r^p^fc cx D{E), 
so that the dependence of the spectrum on the diffusion coefficient cancels out. This implies 
that at the transition between the two regimes one expects a progressive steepening of the 
spectrum of nuclei, asymptotically tending to 6 for a spallation dominated regime. This 
is particularly relevant for heavy nuclei, for which the timescales for spallation and escape 
become comparable in the TeV range. However, the transition region between the two regimes 
is rather broad and in general a low energy hardening of the nuclear spectra (with respect 
to the naive E~^^^) may be expected even if spallation is not dominant over escape. In 
the case of He, with our choice of the parameters {udisc = 3cm~^), Tgp = Tesc occurs at 
Ecrit ~ lOOGe^/n, but the spectrum is sensibly flatter than that of H up to around the 
maximum energy. 

Notice that the disk-like structure of the spatial distribution of sources changes the 
shape of the spectrum of nuclei compared with the simple leaky-box model predictions used 
in Ref. [24]. 

3.1 Contribution of a nearby source 

Here we discuss the relative contribution of an individual source with respect to the diffuse 
CR spectrum in the Galaxy. For simplicity we limit ourselves to the case of protons, while 
the generalization to nuclei of arbitrary charge Z is straightforward. The spectrum of the 
source is assumd to be \]yN{E) where N{E) is the average spectrum injected by SNRs, and 
appearing in Eq. 3.3. The CR density contributed by the (bursting) source is 



{A7rXDD{E)Ty 



where the coefficients \n and \o have been introduced to allow for the possibility that the 
luminosity of the source be somewhat different from the mean luminosity of the CR sources, 
and that the local diffusion coefficient may be different from that averaged over the entire 
volume of the Galaxy. 

This density needs to be compared with the diffuse density as given by Eq. 3.3. The 
condition that the flux from the source dominates over the diffuse CR flux leads to the 
requirement: 

> TT^T^- (3.11) 



{AT,\DD{E)Tf^ - 2^RlD{Ey 
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With this prescription, the condition that a local source dominates the CR flux reads 




( 



R 



) 



5/3 



3GV 



years , 



(3.12) 



where we assumed 7^=1/(30 yr) and Rd = 15 kpc. This result holds within a distance 



from the source. Within such distance the density of CR stays roughly constant. 

The two A parameters have been introduced here with a well defined purpose: the 
luminosity in the form of CRs of a nearby source may be somewhat different from the average 
one, so that A^v allows to consider the possibility that the local source is less powerful {Xn < 1) 
or more powerful {Xn > 1) than average. For simplicity we kept the spectrum of accelerated 
particles the same in shape while only the normalization is allowed to change. 

The meaning oi Xd is more interesting: the diffusion of CRs in the source vicinity can 
be modified by the propagating particles as a result of the self-generation of turbulence. This 
effect leads to lowering the diffusion coefficient with respect to average, thereby leading to 
A^) < 1 (again here we assume for simplicity that only the absolute normalization of D{E) 
changes while retaining the energy dependence). This effect is even more plausible since, as 
discussed above, the number density of CRs within a distance r from the source is larger 
than on average in the Galaxy for a time r, so that self-generation of waves during this time 
may reduce the diffusivity of the particles trying to leave the source. For instance '\i Xd ~ 0.1 
the source contribution dominates on the diffuse CR flux for ~ 100, 000 years (for H = 3 
kpc). Since the energy dependence of r is very weak, this excess is present for a wide range 
of particle rigidities. The distance r is obviously independent of the value of A/). 

4 Estimate of the fluctuations in the spectrum 

One of the goals of the present paper is to discuss the implications of diffusive CR propagation 
from discrete sources on the spectrum measured at Earth. In this section we discuss some 
basic concepts which help us clarifying the role of fluctuations in the source distribution on 
the spectrum. In order to quantify these effects we follow here the formalism of Ref. [27]. 

Based on this formalism, also adopted by [6], given Af independent sources, each one 
producing Nq particles, the average spectrum these produce at a location r at time t is: 



where P{r',t') is the probability of having a source with space-time coordinates {r* ,t') and 
Q{r,t]r',t') is the Green function for transport of particles from {'r',t') to {r,t). For the 
fluctuations in the measured spectrum one finds [27]: 




(3.13) 



{ncR{f,t)) = NoAf / d''r'dt'P{r',t')g{f,t;f',t') 



(4.1) 



{6ncRif,t)6ncR{r' ,t')) 



Sr"dt"P{r",t")g{r,t;T' 



t")Q{^,t'-y, 



t")-\-r{ncR{m{ncR{^,t')). (4.2) 
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In our case, the probability distribution, as already used in Eq. 3.1, is 

nr,t) = ^, (4.3) 

where, by definition, TZ = J\f/At with — )■ oo the time over which the averages are 
calculated. For simplicity, we carry out the calculation for the simple case in which the 
sources are uniformly distributed in an infinitely thin disc with radius R^, and assume for 
our estimates that the spectrum and its fluctuations are measured at the center of the disc, 
r = 0. Moreover we limit our simple estimate to protons, hence ignoring spallation terms in 
the Green function. With these assumptions, the integral over the spatial coordinates reduces 
to an integral over cylindrical radius. The double integral over space and time diverges unless 
a lower bound is imposed on one of the two coordinates, either a minimum distance from 
the closest source, Rmim or a minimum time for receiving its contribution. Train 

: the two are 

obviously correlated, with Rmin ~ V ^DTmin and we choose to impose a Tmin (see [28] for 
a more extended discussion of this effect). The correlation function evaluated in the same 
spatial location is then: 

{6ncnir,t)6ncnir,t)) ^ 3,^3^g,,(^)2r_ ' ^'-'^ 
where we assumed that H'^ / {D{E)Tmin) ^ 1- The strength of the fluctuations is therefore: 

This result shows the strong dependence of 6spec (in terms of normalization and in 
terms of energy dependence) on the cutoff time Tmin, which we imposed in order to avoid 
the divergence in the integral over r. One can see that if T^in is naively chosen as a given 
number, then 6spec turns out to be independent of energy. On the other hand, from the 
physical point of view, the time T^in should carry information about the closest, most recent 
sources around the observer. For instance the time Tmin could be interpreted as the time 
over which one source goes off within a distance from Earth such that CRs from that source 
reach us within a time Tmin- This condition is expressed as 



:!^^[47rD{E)Tm^n] = l^T„ 



AnD{E) 



-1/2 

(4.6) 



With this prescription one obtains [6]: 



Since the diffusion coefficient is chosen as in Eq. 3.4, the expression above can be rewritten 
as: 

Vc(i?) - o.o3F,;f n-,'^'R]%,^^ j , (4.8) 

where we have used D2s/Hkpc ~ 1 (the exact value being irrelevant since the dependence of 
6spec on this ratio is very weak). 
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The energy dependence of the spectral fluctuations is very weak for the standard values 
of 6. At energies E ~ 10^ GeV one can expect fluctuations at the level of ~ 5% for 5 = 1/3 
and at the level of 10% for 5 = 0.6. In practice, as we show below, the fluctuations may be 
somewhat larger because of the combined effect of all nuclei and the occasional dominance 
of some local source contributing most of the CR flux at Earth for some time (see discussion 
in §3.1). This latter scenario is most likely to happen for 5 ~ 0.6. 



5 Results for realistic distributions of SNRs 



In this section we describe our more realistic calculations of spectra and chemical composition, 
taking into account the discreteness in space and time of sources in the Galaxy. First of all 
let us describe our model of CR diffusion in the Galaxy. We assume that the diffusion region 
has a height H above and below the disc, with the latter assumed to have half-width h. In 
the radial direction the diffusion region is not bounded. As we discuss below, this assumption 
is not a serious limitation thanks to the form of the SNR spatial distribution. The diffusion 
coefficient in the whole region is assumed to be spatially constant (the same assumption as 
in propagation codes such as GALPROP) and to depend on rigidity alone, as in Eq. 3.4. 
As suggested by the results of GALPROP, as well as of other propagation codes, a good fit 
to the whole set of available CR data is obtained by taking -D28 /-f^fcpc — 1 (more precisely 
D28/Hkpc = 1-33 for 5 = 1/3 and -D28/^A:pc = 0.55 for 5 = 0.6), therefore we adopt this 
relative normalization of diffusion coefficient and size of the halo. 

We adopt two different models of source distributions. The first one (cylindrical model 
hereafter) simply accounts for the average radial and azimuthal distribution of sources in the 
Galaxy, as discussed for instance in [29]. The other model is instead an attempt at taking 
into account the spiral distribution of sources (hereafter spiral model), for which we adopt 
the formalism of [30]. 

The sources are assumed to be distributed in the Galaxy according to the following 
function [29]: 



■0 



(5.1) 



where /? = 3.53 for Rq = 8.5 kpc. 

The constant A is determined from the normalization condition: 



' dr2TTrf{r) = l^ A=^^^^-^. (5.2) 



The SNR distribution in the z direction is assumed to be 

(5.3) 



f{z) = — exp 



where = 1 is again derived from the normalization condition. 

In the cylindrical model, the positions of the sources are chosen by drawing at random 
values of r and z from the distributions above. The x and y coordinates are chosen by 
generating a random angle {) < cj) < 2ti and using the given value of r. Supernovae are 
generated at a rate TZ = (30 — 100 yr)~^, and the generation of new sources is continued 
until a time span much larger than H^/D[E) is covered, in order to make sure that the 
stationary solution has been reached at the lowest rigidities of interest for us (~ 1 TV). 
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Table 1. Parameters of Galactic arms. 



arm number /name 


K{rad) 


ro(A;pc) 


9o{rad) 


1: Norma 


4.25 


3.48 





2: Carina - Sagittarius 


4.25 


3.48 


vr 


3: Perseus 


4.89 


4.90 


2.52 


4: Crux - Scutum 


4.89 


4.90 


-0.62 



In the spiral model, the procedure we adopt is similar to that introduced in [30]. The 
generation of the position in the z direction is the same as above, therefore we will not discuss 
it any further. A radial coordinate f is drawn at random from the average distribution in 
Eq. 5.1; at this point a random natural number is chosen between 1 and 4, with a flat 
distribution. This number identifies the arm in which the supernova is localized (Norma, 
Carina-Sagittarius, Perseus, Crux-Scutum, as in Table 1). At this point an angular position 
along the arm is associated to the SNR, following the relation: 

9{r) = K log (j^^ +00. (5.4) 

The parameters K, tq and a.re reported in Table 1 (notice that the values of are different 
from those in Table 2 of [30], simply because the axes are rotated by 7r/2 with respect to 
their choice). The Sun is located at {x,y,z) = {8.5kpc,0,0). 

Following the prescription of Ref. [30], we blur the angle 9{r) by Ocorr exp(— 0.35f//cpc), 
where Ocorr is chosen from a flat random distribution between and 27r. Similarly the radial 
coordinate is also blurred by choosing a new value from a normal random distribution with 
mean r and variance 0.07f. A pictorial illustration of the two scenarios is shown in Fig. 1 
where we show the distribution of ~ 30, 000 SNRs in the cylindrical model (left) and the 
spiral model (right). The position of the Sun is illustrated by the thick (red) symbol. 

For each source the spectrum of CRs (protons. He nuclei, CNO nuclei, Mg-Al-Si nuclei 
and Fe nuclei) at the Earth is calculated using the appropriate Green functions, as described in 
§ 2. For each realization of the source distribution in the Galaxy we also compute the chemical 
composition of CRs, as derived from the superposition of the flux of different chemicals. The 
efficiencies of acceleration of nuclei are calculated a posteriori from requiring a fit to the 
available spectra in the TeV region. The calculations are carried out for different choices of 
the propagation parameters. We account for spallation as discussed in § 2, with the average 
gas density in the propagation volume taken as rigas = ridiscih/ H). Notice that with this set 
of prescriptions the grammage traversed by CRs with rigidity R reads 

x(i?)=n,,.cm,^=n,..,/.cm,-^..Xo(3-^j(^^J gem , (5.5) 

which only depends on the density of gas in the Galactic disc and not on the spatial extent 
of the halo, once H^pd ~ 1 is used. In our calculations we adopt UiUsc = 3 cm~^, which 
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Figure 1. A face on view of tlie spatial distribution of SNRs in the Galaxy in the two models: 
cylindrical in the left panel and spiral on the right. About 3 x 10^ sources are shown in each panel. 
Units arc in kpc and the position of the Sun is marked by a thick (red) symbol. 



leads to xo — 9 g cm~^ for 5 = 1/3 and xq ~ 4.3 g cm~^ for 5 = 0.6. It is worth stressing 
that although Eq. 5.5 is normalized at rigidity 3 GV, all of our calculations refer to much 
higher energies, in the TeV range. 

5.1 Spectra and chemical composition for the cyhndrical model 

In this section we discuss many results related to the dependence of the spectra and chemical 
composition on the diffusion coefficient, on discreteness of the sources in space and time and 
the size of the halo for the cylindrical model of source distribution. 

In Fig. 2 we illustrate the all-particle spectrum obtained in 10 realizations of source 
distributions in the cylindrical model, using 5 = \/2>, H = A kpc and a rate of supernova 
explosions in the Galaxy TZ = l/100yr~^ (left) and TZ = l/30yr~^ (right). In both cases 
we impose that the slope 7 of the injection spectrum is related to 5 through 7 + 5 = 2.67, 
where 2.67 turns out to be the slope of the data provided in Ref. [31] as average among all 
experiments (dots with error bars in Fig. 2). The red, staircase line, represents the average 
spectrum resulting from the 10 random realizations. 

As expected, the random distribution of SNRs in space and time leads to fluctuations 
in the all-particle spectrum. One aspect that can be immediately noticed is that the strength 
of the fluctuations is sensibly reduced in the case with a higher SN rate {TZ = l/30yr~^), as 
could be expected based on the simple estimate in Eq. 4.8. Though relatively small, the effect 
of fluctuations is sufficient to induce a visible wiggle in the global shape of the all-particle 
spectrum. 

All spectra present a pronounced knee which compares well with the data and is the 
consequence of the rigidity dependent nature of the acceleration in the sources. In order 
to provide a reasonably good fit to the data one has to assume a maximum energy in the 
proton spectrum Emax = 6 x 10^ GeV, not far from that obtained in non-linear theories of 
particle acceleration at SNR shocks [32]. This number can however reasonably be smaller in 
more realistic models since the spectrum of escaping particles in general does not have an 
exponential cutoff, due to the superposition of advected particles and particles escaping from 
upstream at all times throughout the history of the remnant [33] . The fact that our model 
spectrum underpredicts the data at the highest energies is most likely related to the effect of 
extragalactic cosmic rays, as we discuss below. 
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Figure 2. All-particle CR spectrum for ten random realizations of sources in the cylindrical model, 
assuming 5 = 1/3, a SN rate R = l/lOOt/r^^ [R = l/'dQyr^^) on the left (right). The halo size is 
H = A kpc. The injection spectrum of all chemicals is assumed to have a slope (below the cutoff) 
such that 7 + i5 = 2.67. 

The fluctuations are much more pronounced in the case of faster diffusion (6 = 0.6), as 
illustrated in Fig. 3, an effect that could be expected based on Eq. 4.8. Independent of the 
shape of the average spectrum, one can clearly see that the spectra of individual realizations 
are in general rather wiggly and can fit the observations only in a rather loose sense. The 
bumpiness of the spectrum is the result of the importance of nearby and recent supernova 
events. These realizations also lead to exceedingly large anisotropy as we discuss in Paper 2. 
Both aspects, however, bumpiness and anisotropy, are likely to hint to the breaking of the 
diffusion approximation for CRs coming from nearby sources. 

Using Eq. 3.4 for the diffusion coefficient and recalling that D{E) can also be written 
in terms of the scattering length \{R) as D{R) = (l/3)cA(i?), one easily obtains 

A(i?) p. lOi^Ffcp, cm, (5.6) 

where for simplicity of presentation we assumed D2s/Hkpc = !> independent of 5. 

The diffusion approximation for an individual source is valid only if the source is located 
at a distance d > A, but at R = 10^ GV and with H = A kpc, A w 380 pc for 5 = 0.6. In 
other words at high energy one should be very careful in treating the propagation from nearby 
sources in order to avoid incorrect or artificial conclusions on the spectrum and anisotropy. 

It is interesting to discuss the spectra of individual elements for the different realizations. 
We illustrate our results for 5 = 1/3 for four realizations in Fig. 4. The different lines 
refer to protons (solid/orange). He (dashed/blue), CNO (dot-dashed/cyan), MgAlSi (dor- 
dot-dot-dashed/yellow) and Fe (dotted/green). The red solid line is the all-particle spectrum 
compared with data. 

For 5 = 1/3 the fluctuations are rather small and a good fit to the all-particle spectrum 
can always be obtained. The most important feature shown by these curves is the flattening 
of the spectra of chemicals at lower energies, which results from spallation. In particular the 
spectrum of He appears to be flatter than the spectrum of protons up to ~ 10^ GeV. This 
prediction is observationally consistent with recent CREAM [17, 22] and PAMELA data, and 
implies that at the knee He dominates over protons, even if by a small amount. It is very 
important to realize that the spectral flattening does not imply that spallation dominates over 
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Figure 3. All-particle CR spectrum for ten random realizations of sources in the cylindrical model, 
assuming S = 0.6, a SN rate R = l/30yr~^ and a halo with size H = A kpc. The injection spectrum 
of all chemicals is assumed to have a slope (below the cutoff) such that 7 + (5 = 2.67. 

escape from the Galaxy. The latter condition occurs at lower energies when the escape times 
become larger. For the energies we are interested in, the flattening is simply induced by the 
secular action of spallation during the propagation time, and is certainly more pronounced 
for elements heavier than He, as Fig. 4 clearly shows. In our calculations the slope of the 
proton spectrum is ~ 2.67 but that of the He spectrum is ~ 2.6 (the difference in slope is 
~ 0.05 — 0.07 depending on realization). As shown in Eq. 3.9, if spallation were dominant 
the change of slope would he 5 ^ 0.33. These numbers should be compared with the slopes 
of the proton and He spectra as observed by CREAM, that are 2.66 it 0.02 and 2.58 it 0.02 
[17] respectively. 

The case 5 = 0.6 is somewhat different, as shown in Fig. 5. As stressed above, the fit to 
the all-particle spectrum is not that solid and depends on the specific realization. Moreover, 
the hardening in the He spectrum is not evident and in fact the He and proton spectra are 
mostly parallel in the realizations shown here. Nearby sources induce wiggles and occasionally 
give rise to concave spectra. These are not to be confused with the ones predicted by non 
linear diffusive shock acceleration. 

It is important to keep in mind that the spallation-based explanation of the flatter 
He spectrum compared with the H spectrum relies on the requirement of a large enough 
grammage at energies above 200 GeV/nucleon or so. At these energies the B/C ratio is 
not very constraining. On the other hand, trying to obtain a comprehensive fit of B/C and 
particle spectra at lower energies is also problematic, since the spectra of protons and He 
both show a sharp hardening at 200 GeV/nucleon, that can, equally well, be the result of 
local sources or be due to a reduction of the scattering rate of particles (larger diffusion 
coefficient). The first explanation would not affect the global B/C ratio, while the latter 
would indeed affect it, since it results in a decreased grammage at low energies. In this sense 
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Figure 4. Spectra of chemicals in the cylindrical model, for S = 1/3, a SN rate TZ = l/30yr~^ and 
a halo with size H = A kpc. The injection spectrum of all chemicals is assumed to have slope (below 
the cutoff) such that 7 + (5 = 2.67. The different curves are: solid (orange) for protons, dashed (blue) 
for He, dot-dashed (cyan) for CNO, dot-dot-dot-dashed (yellow) for Mg-Al-Si, and, finally, dotted 
(green) for Fe. Each of the four different panels shows the results for a different random realization 
of the sources' distribution. 

the origin of the harder He spectrum should be investigated together with the origin of the 
change of spectral slope seen for all species at 200 GeV/nucleon. 

From the point of view of energetics the two scenarios with (5 = 1/3 and 5 = 0.6 are 
appreciably different. On average the case 5=1/3 leads to require an efficiency of acceleration 
(for protons) of the order of ~ 6 — 7% for TZ = l/30yr~^. For 5 = 0.6 the efficiency of proton 
acceleration must be ~ 10%. The fact that (especially for 5 = 1/3) the acceleration efficiency 
for individual SNRs is so small might be a possible reason why it is very difficult to detect 
the gamma ray emission from pion production and decay. It can also provide a useful hint 
to the development of non linear theories of shock acceleration: on one hand the acceleration 
efficiency is not required to be high, on the other it has to be large enough to induce the 
magnetic field amplification necessary for particle acceleration up to the knee region. 

A summary of the results on chemical composition is provided by the mean value of the 
logarithmic mass as a function of energy. We plot this function in Fig. 6 for 5 = 1/3 (left) 
and 5 = 0.6 (right). The results of our calculations are compared with the data compiled in 
Ref. [24]. Unfortunately the large error bars in the data do not allow to add information to 
what we have already pointed out above. The behaviour of the data in the highest energy bins 
might be interpreted as suggestive of the need for a lighter component, possibly contributed 
by extragalactic cosmic rays. 
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Figure 5. Same as for Fig. 4 but for 6 = 0.6. 




Figure 6. Mean Log of the mass of CRs for ten random realizations of the distribution of SNRs in 
the cyhndrical model, assuming S = 1/3 {S = 0.6) on the left (right) and a halo with size H ^ 4 kpc. 
The red staircase line is the average over the random realizations. 

5.2 Spectra and chemical composition for the spiral model 

In this section we discuss the implications of the spiral distribution of SNRs in the Galaxy 
for the spectrum and chemical composition of cosmic rays observed at Earth. As shown in 
Fig. 7, again a very good fit to the measured all-particle spectrum can be obtained adopting 
the same conditions as for the cylindrical model, namely 7 + (5 = 2.67. In this case, however, 
an interesting dependence of the results on the size of the magnetized halo arises. 

From Eq. 5.5 we see that the grammage depends very weakly on H because of the 
condition Hkpc/D2s ~ 1- This implies that changing (decreasing or increasing) H has no 
appreciable consequences in the case of the cylindrical model, where the source density around 
the Sun is roughly constant (within a distance of order ~ H). This conclusion remains 
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Figure 7. Spectra of chemicals in tlie spiral model, for S ~ 1/3, a. SN rate TZ = l/30yr~^ and a halo 
with size H = 4 kpc. The injection spectrum of all chemicals is assumed to have slope (below the 
cutoff) such that 7 + i5 = 2.67. 

approximately true also in the spiral model, although here changing H leads to a change 
in the volume (and number of arms) contributing a flux at Earth. Now, by taking a look 
at Fig. 1, it is immediate to notice that the position of the Sun is right in between two 
Galactic arms, so that most of the sources are some distance away. Nuclei that have to be 
transported from the arm to the Sun suffer spallation and their depletion at low energy is not 
compensated by the contribution of nearby sources. This is exactly the situation in which 
changing the volume of contributing sources might make a difference. 

This fact is illustrated in Fig. 8 where we show the all-particle spectrum for H = 2 
(left) and H = 4 kpc (right). For H = A kpc, we have already seen that the spiral model 
reproduces the data with a spectral index at injection such that 70^^ = 2.67: the same result 
that we obtained for the cylindrical model, independent of H. On the other hand, for H = 2 
kpc, 7o6s = 2.67 is only marginally compatible with the data, while 7 + 5 = 2.7 is found to 
be favorite. This is perfectly consistent with the fact that in the spiral model, for a given 
value of 7 + (5, heavy nuclei have a slightly harder spectrum. The difference in slope between 
the two cases, even if small, makes a sensible difference on the all-particle spectra because it 
sums up over many decades in energy. 

In perfect analogy with the cylindrical case, the secular action of spallation leads also 
in this case to the flattening of the spectra of heavy elements for 5 = 1/3. The He spectrum 
is again flatter than the protons' by the same amount as in the previous model, and again 
up to an energy of ~ 10^ GeV. 

5.3 Comments on the He hardening 

Two important observational findings recently appeared in the literature on cosmic rays: 
first, the sharp hardening found by ATIC-2 [34], CREAM [17] and PAMELA [35] in both the 
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Figure 8. Dependence of the all-particle spectrum on the realization of sources in the spiral model, 
for H = 2 (left) and H = A (right) kpc. In the left (right) panel we used 7 + (5 = 2.7 (7 + 5 = 2.67). 



proton and He spectra at ~ 200 GeV/nucleon; second, the fact that the He spectrum above 
this energy appears to be harder than the proton spectrum. The latter effect is observed at 
a very high level of significance by all three experiments mentioned above. 

Taken at face value, the first finding alone represents a serious challenge to our attempts 
at achieving a comprehensive fit to all cosmic ray observables, an attempt that is routinely 
made by using GALPROP or similar propagation codes. Two possibilities come to mind: that 
the steep, low energy, part of the spectrum may be contributed by a so-far unknown nearby, 
probably old, cosmic ray source; or that at low energy cosmic rays experience a scattering rate 
that decreases with decreasing energy more than expected based on extrapolation from high 
energies, namely the diffusion coefficient has a stronger energy dependence with consequent 
steepening of the equilibrium spectra. 

The first scenario effectively decouples the B/C ratio from the measured spectral shape 
of primary cosmic rays that produce B through spallation reactions: the production of B 
depends on the mean Galactic CR spectrum, which in this case would be different, at low 
energies, from the spectrum at Earth. In this scenario it is virtually impossible to know 
the actual spectrum of CRs in the Galaxy, though some clues can be gathered from diffuse 
backgrounds, such as gamma rays and radio emission. It is worth stressing that in this case, 
propagation codes such as GALPROP cannot be easily used since they do not include the 
possibility of taking into account discrete sources. 

In the second scenario one expects the dependence on energy of the grammage to be 
stronger at low energies than it is at high energies, thereby leading to a steepening of the 
equilibrium spectra of all cosmic ray components. The B/C ratio would reflect this complex 
situation, but would not provide any strong constraint on the grammage experienced by high 
energy particles. 

In both situations extrapolating functional forms from high to low energies appears 
unjustified and potentially dangerous in that it may lead to flawed conclusions. In a recent 
paper [36], among other things, the authors also comment on the results presented on a 
preliminary version of the present paper (which had been made publicly available as an 
arXiv preprint): the authors extrapolate our curves and conclusions to very low energies 
in order to compare predictions with B/C data where they are available. They conclude 
that if one interprets the He harder spectrum as a consequence of spallation, the grammage 
traversed by particles at lower energies leads to exceed the fiux of antiprotons and the B/C 
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ratio. However, antiproton measurements are only available at energies lower than 200 GeV 
[37], namely outside the range of applicability of our model, which at those energies cannot 
even reproduce the CR spectrum, as already stressed above. On the other hand, for the B/C 
ratio there are two measurements, by CREAM [38] and TRACER [39] at energies higher 
than 1 TeV, but both have large error bars. If these two points are taken with the error bars 
shown in [39], no inconsistency can be claimed between our model and the relevant data. 
The obvious conclusion is then the same that could be reached by observation of the spectral 
steepening suffered by all chemicals below ~200 GeV/nucleon: the grammage cannot be 
extrapolated down in energy with the same slope that we assume at TeV energies. 

As far as the spallation-related explanation of the He hardening is concerned, we think 
that there is still enough uncertainty in the B/C ratio in the TeV range to warrant the 
benefit of doubt on whether this may be the correct explanation. In our view this scenario 
is certainly less ad hoc and more physically motivated than assuming artificial breaks in the 
source spectra or unknown source populations. To our knowledge this is the only attempt 
so far at understanding the new data on differential hardening in terms of well established 
propagation physics, rather than simply fitting them. 

Finally, a very important feature of this scenario is that it has the potential to be 
disproved not too far in future, through more accurate measurements of secondaries in the 
high energy range, with TRACER and AMS-2. While extrapolation to energy regions where 
the initial assumptions are known not to hold might lead to misleading conclusions, we think 
that any attempt should be made to test whether the spallation of He at TeV energies may 
lead to observables that are in contradiction with data in the same energy range. 

6 Additional scenarios 

Other physically relevant scenarios have been tested but did not induce particularly interest- 
ing effects on the spectrum and chemical composition: one is the continuous leakage of CRs 
from SNRs, as discussed in § 2, and the other is the scenario in which supernova explosions 
occur dominantly in OB associations. Both these assumptions have negligible consequences 
in terms of changes induced in the spectrum and chemical composition of CRs observed at 
Earth. 

7 The effect of extragalactic cosmic rays 

The spectrum and the chemical composition in the energy region above the knee may be 
substantially affected by a contribution from CRs of extragalactic origin. This is true in both 
descriptions of the transition from galactic to extragalactic that are presently considered as 
most likely, namely the electron-positron dip model [40] and the mixed composition scenario 
[41]. Since the investigation of the transition is not among the purposes of the present paper, 
here, for simplicity, we illustrate this effect only for the case of the dip. We calculate the 
spectrum of extragalactic CRs as discussed in detail in [42] and normalize the absolute flux 
to the HiRes data ([43] and references therein). The injection spectrum of extragalactic CRs 
is assumed to be a power law with slope —2.7. No evolution of the source luminosity with 
redshift is taken into account since the only purpose of the calculation is to illustrate the 
effect, and not to explore the region of allowed parameters. As discussed in [42], the dip 
model provides a suitable explanation of the observations provided that no more than ~ 15% 
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Figure 9. All particle spectrum of CRs including the extragalactic component with Ecut = 10^ (left) 
and Ecut = 10*^ GeV (right). 

of the extragalactic flux is contributed by He nuclei. Therefore we assume a pure proton 
composition of extragalactic CRs. 

Even a small magnetic field in the intergalactic medium can suppress the flux of CRs 
coming from extragalactic sources, thereby introducing a low energy cutoff in the spectrum. 
Such low energy cutoff can also be intrinsic to the source if the sources harbor relativistic 
shocks: in this case the minimum energy of the accelerated particles is ~ AT'^mpC^, with T 
the Lorentz factor of the shock. 

We artificially model this low energy suppression with an exponential cutoff at Ecut = 
lO'^ GeV and 10^ GeV. The all-particle spectrum for these two values of E^ut is plotted in 
Fig. 9: one can see that the superposition of the Galactic and extragalactic CR spectra 
allows us to provide a smooth fit to the available data in both cases. In Fig. 10 we plot 
the mean logarithmic mass, (log A), as a function of energy for the two values of Ecut- the 
goodness of the fit is negatively affected by an extragalactic proton spectrum extending too 
low in energy, therefore if we are to believe this type of interpretation of the transition, and 
take strictly the data on the mean log mass, then the low energy cutoff of the extragalactic 
CR spectrum has to be at > 10*^ GeV. At very high energy, Fig. 10 suggests that the 
transition has to occur through mixing of different chemicals. However one has to keep in 
mind that the measurements of the composition in the transition region are still affected 
by large systemic errors which do not allow very firm conclusions in this important energy 
range, as also discussed in [44]. Recent data from HiRes [43] and Telescope Array [45] seem 
to find a chemical composition at ~ 10^^ eV dominated by protons, which is not reflected 
in the data used in Fig. 10. The most recent data on the chemical composition as observed 
with the Pierre Auger Observatory also show a proton dominated composition at ~ 10^*^ eV, 
but gradually shifting to heavier composition at higher energies (see for instance [46]). 

Interestingly enough our results do not require an additional component (possibly super- 
heavy Galactic nuclei invoked by [47] and [31], see also [48]) to fit the shape of the all-particle 
spectrum. The need for such an extra component was probably arising from assuming that 
the maximum energy of protons accelerated in SNRs equal the energy where the knee appears. 
In our calculations, as discussed above, the maximum energy of protons that best fits the all 
particle spectrum is ~ 6 x 10^ GeV, well above the knee, but the knee itself appears to be 
mainly contributed by Helium. The shape of the cutoff in the spectra of individual chemicals 
also affects the need for extra-components. 
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Figure 10. Mean value of log A including the extragalactic component with Ecut = lO*" (left) and 
E,ut = 10« GeV (right). 

8 Summary and Conclusions 

In this paper we have tried to make progress in the quest for the origin of CRs, trying to 
clarify the requirements on propagation and source spectra that come from observations of 
CRs at Earth. 

We computed the spectrum and chemical composition of CRs at Earth assuming that 
SNRs are the primary sources of galactic CRs and taking into account the discrete distribution 
in both space and time of these sources. Each source is assumed to inject protons and 
nuclei with a power law spectrum up to some maximum energy, above which the spectrum is 
exponentially cut off. The power law index and maximum energy are assumed to be the same 
for all sources and are determined from observations, with Emax simply scaling with rigidity 
for nuclei heavier than H. The propagation of CRs from their sources to Earth is treated 
through a Green function approach allowing for a completely general dependence on space 
and time of the injected spectrum, as well as taking into account spallation. The diffusion 
coefficient has been taken constant in space within the galactic halo and has a power-law 
dependence on rigidity alone with index 6 = 1/3 or S = 0.6. In terms of its normalization, 
the diffusion coefficient is assumed to scale with the size of the halo, following the findings 
of other propagation codes, such as GALPROP or DRAGON. 

We have considered a few different scenarios for the distribution of sources. As far as 
their spatial distribution is concerned, we have carried out the calculations both for a smooth 
distribution of SNRs within the galactic disk (cylindrical model) and for the case in which 
SNRs follow the spiral structure of the Galaxy (spiral model). We have also considered the 
possibility that SNRs are clustered in superbubbles but we find that in terms of spectrum 
and chemical composition this assumption makes negligible difference. 

Concerning the dependence on time of the injected spectrum, we have considered two 
different scenarios: a case in which all the CRs a SNR has ever accelerated are released at 
the end of its evolution, and a case in which, at all times through the Sedov phase of the SNR 
evolution, there is continuous leakage of CRs at the maximum energy from upstream of the 
decelerating blast wave, with the maximum energy decreasing with time down to ~ IGeV. 
These two different models for the time dependence of the injection lead to very similar 
outcomes in terms of spectrum and chemical composition, once the spatial distribution of 
the sources and the energy dependence of the diffusion coefficient are fixed. 



Our main findings in the present paper are summarized in the following. First of all 
the ability of our approach at taking into account the discreteness of the source distribution 
allows us to illustrate that the all-particle CR spectrum at Earth is not in general guaranteed 
to be a power law. In particular if the diffusion coefficient depends on rigidity as with 
S = 0.6 a number of wiggles are likely to appear in the spectrum, due to the contribution 
of nearby sources. This fact also suggests that this value of 5 will lead to large levels of 
anisotropy. 

For a diffusion coefficient scaling with E^^^, the CR spectrum and chemical composition 
observed at Earth can be reproduced in a satisfactory way in both the cylindrical and spiral 
model of source distribution in the Galaxy. The required maximum energy of the accelerated 
particles is Emax = 6 x 10^ GeV for protons, reasonably close to (but somewhat higher 
than) what can be expected from a SNR based on the non linear theory of diffusive shock 
acceleration. The knee in the all-particle spectrum results from the scaling of the maximum 
particle energy with rigidity. The required power-law index of the injected particles is the 
same in both scenarios for a halo size H = 4 kpc, namely 7 = 2.34. A slightly steeper 
power-law is required instead in the spiral case if the magnetized halo has size H = 2 kpc. 
In this latter case one finds 7 = 2.37, and the small difference is induced by the need of 
compensating for the effects of spallation that lead to a hardening of the spectrum at low 
energy, with a deficiency of low energy particles that in this scenario cannot be supplied by 
nearby sources due to the inter-arm position of the solar system. 

These spectral indices represent at present the most problematic aspect of the supernova 
remnant paradigm for the origin of CRs. The required acceleration efficiency is relatively 
low, for D cx E^/^ and TZ = 1/30 yr~^ being only of order 6 — 7%. 

Proper account of the effects of spallation leads to a hardening of the spectrum of heavy 
elements at low energy, even when this process is not dominant over escape. This is the key 
element that leads us to expect that, for 5 = 1/3, He nuclei have a flatter spectrum than 
protons at low energies and dominate the all-particle spectrum at the knee, in qualitative 
agreement with the recent findings of CREAM [17, 22] and PAMELA [35]. The slopes 
of the protons and Helium spectra that we find are ~ 2.67 and ~ 2.6 respectively, to be 
compared with the measured values of 2.66 it 0.02 and 2.58 it 0.02 [17]. Our calculations show 
that for 5 = 0.6 this effect disappears, that might indeed suggest that Galactic diffusion is 
characterized by a low value of 6. 

The prediction of the hardening in the spectra of nuclei compared with that of protons 
is one of the most important results of our calculations. As noted by [36], if extrapolated 
down to much lower energies this finding may cause excessive B and antiprotons production. 
However, a more straightforward consequence of such a naive extrapolation would be that 
the spectra of H and He would not reflect the observed steepening of both components at 
rigidities < 200 GeV/nucleon, which we make no attempt at explaining. In other words, 
what this is suggesting is simply that at low energies something new is happening. Currently 
this is not yet explained in physical terms: the data can only be reproduced by introducing 
artificial breaks in either the injection spectra or the diffusion coefficient, or both. We limited 
ourselves to consider the high energy range and the simplest possible physical scenario for the 
very reason that we did not want to introduce additional parameters before having an idea of 
the nature of the underlying phenomena. After all, if one applied the same approach to the 
knee in the CR spectrum, assuming breaks in the injection spectra and/or in the diffusion 
coefficient could certainly lead to a good fit to the data, but would make little contribution 
to our understanding of the origin of this spectral feature. 
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Currently we cannot explain the fact that the spectra of all species flatten rather 
abruptly at ~ 200 GeV/n as found by the PAMELA experiment [35]. However we do not 
share the opinion of the authors claiming the failure of the SNR paradigm or the existence 
of an unknown accelerator. We rather think that the SNR paradigm is in fact more complex 
than usually assumed in doing these claims, and that its consequences are not yet so well 
understood as sometimes people would like to believe. 

At high energies, above the knee, we underpredict the all-particle spectrum and slightly 
overpredict the mass mean Log of CRs. This seems to suggest that a lighter CR component is 
probably contributed at those energies by extragalactic sources. Since fitting the details of the 
observed CR spectrum at high energy was not our main purpose in this paper, we modeled this 
extragalactic component only within the scenario of the Dip [40], in which the extragalactic 
CR composition is dominated by protons. In this case we find that the spectrum can be fitted 
very well, while the chemical composition suggests that the extragalactic component should 
have a low energy cutoff at energies larger than 10^ GeV. Fitting the all-particle spectrum 
does not require any additional heavy component as previously postulated in Refs. [47] and 
[31]. 

Our general conclusion is that if the diffusion coefficient in the Galaxy scales with 
energy as £■1/3 SNRs can account for the all-particle spectrum and chemical composition of 
CRs detected at Earth, provided that they are able to accelerate protons up to ~ 6 x 10^ 
GeV and together give spectra as steep as E~'^ with 7 = 2.3 — 2.4. This is indeed a bit 
of a challenge for the non-linear theory of diffusive shock acceleration, which is believed to 
describe the acceleration process in this context. On the one hand, efficient acceleration 
leads to hard (and even concave) particle spectra, at odds with the relatively large values of 
7 obtained above. On the other hand, the inferred maximum energy can only be achieved if 
the magnetic field is largely amplified with respect to its average value in the ISM. Efficient 
amplification is usually associated with efficient particle acceleration, whereby the difficulties 
at reconciling steep spectra and high maximum energies. These difficulties might be overcome 
in the context of NLDSA if the velocity of the scattering centers that enters the transport 
equation is that computed in the amplified magnetic field, rather than the unperturbed Alfven 
velocity [49]. In fact, the use of the unperturbed velocity derives from a strict application of 
quasi-linear theory, which is however questionable in the presence of such large levels of field 
amplification as those usually found appropriate for young SNRs. In addition, as soon as the 
Alfven velocity in the self-generated magnetic field is considered, the acceleration efficiency is 
much reduced, while it may still happen that the maximum proton energy reaches the knee. 
This recent result is very interesting, especially in the light of our findings in this paper, that 
the required acceleration efficiency is only about 5% in our best fit scenario, suggesting that 
the acceleration process has to guarantee decoupling between a large maximum energy and 
efficient particle acceleration. 

Another possibility to reconcile steep spectra of the accelerated particles and a high 
maximum energy might be related to the shock obliquity. The theory of NLDSA has so far 
been developed only for parallel shocks, while particle acceleration at perpendicular and more 
generally oblique shocks has not been devoted the same amount of attention. Perpendicular 
shocks are fast accelerators, and in principle allow to reach higher energies than parallel 
shocks. This is especially interesting in that the value we require for the maximum energy 
of protons is close to the upper limit of what can be achieved within NLDSA at parallel 
shocks. Whether the required combination of steep spectra and high maximum energies can 
be achieved at highly oblique shocks is an issue that deserves further investigation. 
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